{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6b95dff3830dd37f",
   "metadata": {},
   "source": [
    "事先说明，以下绘图数据皆为虚拟数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4126a2f4-9634-4177-91e4-c75fc8094ec2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.font_manager import FontProperties\n",
    "import matplotlib as mpl\n",
    "\n",
    "# 设置字体路径\n",
    "# linux 下\n",
    "# font_path = \"/usr/share/fonts/google-noto-cjk/NotoSansCJK-Medium.ttc\"\n",
    "# font_prop = FontProperties(fname=font_path)\n",
    "# # 将字体设置应用到matplotlib的配置参数中\n",
    "# mpl.rcParams['font.family'] = 'sans-serif'\n",
    "# mpl.rcParams['font.sans-serif'] = [font_prop.get_name()]\n",
    "\n",
    "# windows 下\n",
    "font_path = \"C:\\\\Windows\\\\Fonts\\\\msyh.ttc\"\n",
    "font_prop = FontProperties(fname=font_path)\n",
    "plt.rcParams['font.family'] = 'sans-serif'\n",
    "plt.rcParams['font.sans-serif'] = [font_prop.get_name()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94a9238e-0977-4882-9633-5452f4d06d68",
   "metadata": {},
   "source": [
    "# 引言\n",
    "\n",
    "近年来国家气象部门大力发展遥感观测系统，诸如这几年建设的地基垂直观测系统包含：风廓线雷达、微波辐射计、云雷达、激光雷达、GNSS/MET 观测设备，这些设备实现了对高空大气 24 小时不间断的自动观测，解决了高空大气探测数据少的问题。\n",
    "上一章节我们学习风廓线雷达探测原理和数据格式，绘制了风羽图；\n",
    "\n",
    "本次我们一起学习微波辐射计的采集原理，并利用 metpy 绘制 TlogP 图。\n",
    "\n",
    "在学习完新设备的原理、数据处理和绘图后，我们再以 闪电定位观测数据为例，一起学习如何处理大量数据。\n",
    "\n",
    "#  案例 1：利用微波辐射计数据绘制 TLogP 图，并叠加风廓线雷达\n",
    "\n",
    "微波辐射计（ microwave radiometer，缩写为 “MWR”）也称为 “微波辐射仪”，是一种用于测量亚毫米级到厘米级波长（频率约为 1-1000GHz）的电磁波（微波）的辐射计。微波辐射仪能接收大气中的某些成分在一定频率上强烈辐射的微波，经过一定的转换方法，得到大气在垂直和水平方向上的气象要素分布。此仪器携带方便，可增加探空网在时间和空间上的密度，能观测到大气的连续变化，不致漏掉范围较小但变化剧烈的天气系统。\n",
    "\n",
    "最早的微波辐射计是由美国物理学家罗伯特・亨利・迪克于 1946 年发明的，我国微波辐射计的研究工作始于 70 年代初，已经历了从地基到空载，从单频到多频的发展进程，并取得一系列有价值的成果。\n",
    "下图为暴雨所微波辐射计安装现场图 ，图片来源于 湖北气象官网\n",
    "\n",
    "![国内微波辐射计](https://cdn.kesci.com/upload/image/rl0n6yavao.jpg?imageView2/0/w/960/h/960)\n",
    "\n",
    "\n",
    "## 探测原理\n",
    "\n",
    "地基微波辐射计是基于物质热辐射理论进行大气探测，它接收特定频段的辐射，经过数值化处理得到有效的信息，这些频段的辐射电压值就是 L0 级数据；\n",
    "对辐射电压值进行按照普朗克辐射定律进行数学处理，得到 L1 级数据 —— 亮温；\n",
    "再通过人工神经网络算法对 L1 级的亮温数据进行反演，得到 L2 级数据 —— 高空大气 的：气温、湿度、水汽密度、液态水含量等数据。\n",
    "\n",
    "探测数据采集处理流程如下：\n",
    "微波多通道观测→L0 级数据（电压值）→L1 级数据（亮温）→L2 级数据（大气廓线）\n",
    "\n",
    "日常业务数据处理中很少使用到 L0、L1 级数据，一般使用 L2 级大气廓线数据。\n",
    "\n",
    "## 大气廓线数据格式\n",
    "\n",
    "微波辐射计 L2 级大气廓线数据采取 txt 文本 和 bin 二进制两种方式存储，一般业务使用 txt 文本数据格式，该格式比较简单，具体如下：\n",
    "\n",
    "> MWR,01.00\n",
    "57494,114.0600,30.6000,+23.6,TQ9527,lt\n",
    "Record,DateTime,10,SurTem (℃),SurHum (%),SurPre (hPa),Tir (℃),Rain,CloudBase (km),Vint (mm),Lqint (mm),0.000 (km),0.025 (km),0.050 (km),0.075 (km),0.100 (km),0.125 (km),0.150 (km),0.175 (km),0.200 (km),0.225 (km),0.250 (km),0.275 (km),0.300 (km),0.325 (km),0.350 (km),0.375 (km),0.400 (km),0.425 (km),0.450 (km),0.475 (km),0.500 (km),0.550 (km),0.600 (km),0.650 (km),0.700 (km),0.750 (km),0.800 (km),0.850 (km),0.900 (km),0.950 (km),1.000 (km),1.050 (km),1.100 (km),1.150 (km),1.200 (km),1.250 (km),1.300 (km),1.350 (km),1.400 (km),1.450 (km),1.500 (km),1.550 (km),1.600 (km),1.650 (km),1.700 (km),1.750 (km),1.800 (km),1.850 (km),1.900 (km),1.950 (km),2.000 (km),2.250 (km),2.500 (km),2.750 (km),3.000 (km),3.250 (km),3.500 (km),3.750 (km),4.000 (km),4.250 (km),4.500 (km),4.750 (km),5.000 (km),5.250 (km),5.500 (km),5.750 (km),6.000 (km),6.250 (km),6.500 (km),6.750 (km),7.000 (km),7.250 (km),7.500 (km),7.750 (km),8.000 (km),8.250 (km),8.500 (km),8.750 (km),9.000 (km),9.250 (km),9.500 (km),9.750 (km),10.000 (km),QCflag\n",
    "1,2022-10-29 00:00:00,11,13.60,100.00,995.20,13.62,1,0.00,45.15,3.28,13.59,10.62,9.95,10.87,9.43,9.08,8.85,7.13,7.48,9.84,8.23,9.77,10.54,9.52,9.19,8.56,10.75,10.23,11.28,9.91,11.15,10.52,11.74,12.13,12.68,12.31,13.99,16.06,14.42,15.3,16.8,15.74,15.22,17.07,18.16,17.46,16.56,18.54,18.13,19.28,20.81,19.5,19.99,21.35,21.46,21.64,20.94,21.24,20.66,19.85,20.16,20.67,17.29,15.0,13.13,13.23,12.14,10.1,7.87,8.72,5.42,4.88,4.65,3.01,0.04,-0.79,-2.82,-4.5,-7.83,-9.0,-9.14,-9.59,-11.9,-13.43,-14.96,-17.23,-17.89,-20.09,-24.16,-27.42,-31.12,-32.7,-36.28,2\n",
    "2,2022-10-29 00:00:00,12,13.60,100.00,995.20,13.62,1,0.00,45.15,3.28,11.08,16.48,16.9,17.4,18.59,18.28,19.49,19.34,19.38,19.29,19.1,19.06,20.31,19.87,20.57,19.65,20.04,20.74,20.12,21.26,20.27,20.91,20.79,19.87,20.64,20.27,21.01,20.49,19.16,20.18,19.88,19.67,19.62,20.02,19.46,19.67,19.34,19.32,18.41,18.38,18.57,19.05,17.95,17.51,16.81,17.18,16.4,16.88,16.56,15.36,15.43,13.75,10.48,7.21,5.02,3.43,2.13,0.86,0.79,0.71,0.31,0.11,0.0,0.0,0.0,0.17,0.27,0.23,0.25,0.25,0.49,0.37,0.39,0.5,0.48,0.45,0.4,0.38,0.35,0.33,0.28,0.25,0.22,2\n",
    "3,2022-10-29 00:00:00,13,13.60,100.00,995.20,13.62,1,0.00,45.15,3.28,100.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,95.14,86.19,73.01,58.48,51.71,43.62,42.52,43.77,40.18,38.84,36.93,35.78,40.7,34.08,30.75,35.0,30.86,29.9,31.54,32.07,36.88,32.42,29.41,39.57,42.06,40.46,41.22,35.79,45.41,53.18,2\n",
    "4,2022-10-29 00:00:00,14,13.60,100.00,995.20,13.62,1,0.00,45.15,3.28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\n",
    "5,2022-10-29 00:02:00,11,13.60,100.00,995.20,13.75,1,0.00,44.78,3.36,13.59,12.18,11.63,12.73,11.63,11.65,11.59,10.08,10.45,12.63,11.08,12.49,12.89,11.91,11.74,10.99,12.86,12.33,13.01,11.57,12.67,11.89,12.74,12.84,13.2,12.68,13.87,15.28,13.72,14.37,15.59,14.39,13.92,15.33,15.86,15.38,14.31,15.96,15.7,16.57,18.06,16.64,17.16,18.21,18.36,18.57,18.01,18.19,17.76,16.88,17.5,18.18,15.46,13.82,12.29,12.73,11.82,10.04,8.13,8.99,5.77,4.9,4.65,2.99,0.03,-0.97,-3.05,-4.81,-7.83,-9.25,-9.77,-10.33,-12.94,-14.71,-16.25,-18.63,-19.41,-21.5,-25.59,-28.68,-32.04,-33.66,-36.89,2\n",
    "\n",
    "第一行为标识字符串 ，**MWR 表示为 微波辐射计数据** , 01.00 表示为数据版本号为 1.0\n",
    "第二行为站点元数据信息，使用 `,` 分隔 ，该行数据以此为：\n",
    "\n",
    "|  序号 | 名称  | 说明  |\n",
    "| ------------ | ------------ | ------------ |\n",
    "|  1 |  站号 |  ------  |\n",
    "|  2 |   站点经度 | ------  |\n",
    "|  3 |   站点纬度 | ------  |\n",
    "|  4 |   站点拔海高度 | ------  |\n",
    "|  5 |   探测仪器型号 | ------  |\n",
    "|  6 |   型号标识 | ------  |\n",
    "\n",
    "从第三行起到结尾，内容类似 CSV 格式 ，可跳过首两行后，使用 pandas.read_csv 读取，各列意义如下：\n",
    "\n",
    "|  序号 | 列名  | 内容说明  |\n",
    "| ------------ | ------------ | ------------ |\n",
    "|  1 |  Record |  记录序号  |\n",
    "|  2 |   DateTime | 数据观测时间  |\n",
    "|  3 |   10 | 廓线数据类型  |\n",
    "|  4 |  SurTem (℃) | 环境温度（单位为℃）  |\n",
    "|  5 |   SurHum (%)  | 环境相对湿度（单位为 %）  |\n",
    "|  6 |  SurPre (hPa) | 大气压力（单位为 hPa）  |\n",
    "|  7 |  Tir (℃)  |   红外温度（单位为℃）     |\n",
    "|  8 |     Rain  |  是否降水（1＝是，0＝否）    |\n",
    "|  9 |  CloudBase (km)  | 云底高度（单位为 km）|\n",
    "|  10 |     Vint (mm)   |  积分水汽（单位为 cm） |\n",
    "|  11 |    Lqint (mm)    | 积分云液水（单位为 mm） |\n",
    "| 12... |        Hn (km))    |   第 n 层结数据，每层一列，不同设备层结数量不同 |\n",
    "\n",
    "其中 `10` 列 表示 L2 级 气象产品数据组廓线数据的不同类型，具体意义如下：\n",
    "\n",
    "|  类型码  |   廓线数据类型  |   层结数据单位   |\n",
    "| ------------ | ------------ | ------------ |\n",
    "|  11  |   温度廓线 TEM   |  ℃    |\n",
    "|  12   |  水汽密度廓线 WVDen   | $g/m^{3}$   |\n",
    "|  13   |  相对湿度廓线  RHU  | %      |\n",
    "|  14   |  液态水廓线 WDen  | $g/m^{3}$   |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd3be5b-cbd9-4bd4-838c-a47821388c5b",
   "metadata": {},
   "source": [
    "## 与气球探空相比\n",
    "\n",
    "微波辐射计的观测频率为 **2 分钟一次**，每次观测可以得到 **4 种廓线产品**，一天共有 **2880 条数据**。\n",
    "\n",
    "而每日探空站只能观测两次，从观测时间频率上 微波辐射计有着较好的观测密度，能更好地监测灾害性天气的发生与发展。\n",
    "\n",
    "但微波辐射计数据是遥感反演高空数据，使用 人工神经网络算法建立了 亮温 到 大气廓线的映射关系，在反演中存在着误差。\n",
    "\n",
    "目前在实际业务中还是以气象探空为主要高空气象数据来源，微波辐射计等垂直观测系统通过数据融合方式对起到对高空数据补充。\n",
    "\n",
    "（世界气象组织定期组织各成员国对气球探空仪器进行对比实验，2010 年以后我国探空仪器设备研发逐渐追上国际脚步，数据质量逐年提高，特别是在 2018 年气象仪器设备厂商整合后，数据质量已经与国际水平不相上下。）\n",
    "\n",
    "## 数据读取示例"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb8f4972-fd86-4510-9a07-677a279e043a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "sys.path.append(os.path.abspath('../'))\n",
    "import json\n",
    "from datetime import datetime,timedelta\n",
    "import dateutil.parser\n",
    "import dateutil.rrule\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "from metgrs import Utils\n",
    "import io\n",
    "import math\n",
    "import types\n",
    "import glob\n",
    "import matplotlib as mpl\n",
    "from joblib import Parallel,delayed\n",
    "from metgrs import MicroWaveRadiometer,WindProfileRadar\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c290e167-2dac-489a-b9f2-08dec603f4cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath='../datas/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc155483-60a8-4b8c-a981-a91671ce0f83",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs=sorted(glob.glob(datapath+'/MWRFiles/20241204/*'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28966ceb-0bd5-4d70-bfa8-bb8e868b151d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=MicroWaveRadiometer.readMWRFile(fs[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b29ec9f-2881-4989-8284-2d479f13fb7b",
   "metadata": {},
   "source": [
    "## TLogP 图\n",
    "\n",
    "又名 SkewT-LogP 图，是一个非常有用的气象热力学图，在该图上绘制了某个站点高空大气的不同温度、湿度、气压等气象要素，图的横坐标为 气象要素的分布，纵坐标为 气压的对数。\n",
    "\n",
    "\n",
    "![TLogP](https://cdn.kesci.com/upload/image/rl0ruy1fx6.png?imageView2/0/w/960/h/960)\n",
    "\n",
    "TLogP 图上常常会绘制三条曲线：\n",
    "\n",
    "- 温度层结曲线 ，又称为大气温度廓线，表明真实大气不同高度的气温。\n",
    "\n",
    "- 露点层结曲线，又称为大气露点廓线，表明真实大气不同高度的露点温度。\n",
    "\n",
    "- 状态曲线，这是一个虚拟的曲线，用于模拟地面气块沿着干绝热上升达到饱和后，再沿着湿绝热上升的气块温度。\n",
    "\n",
    "三条曲线结合相关气象指标，在气象灾害预报中有着重要的作用。\n",
    "\n",
    "- 是暴雨、冰雹等强对流天气预测的重要指标；\n",
    "\n",
    "- 是云高度和厚度的重要识别方法\n",
    "\n",
    "- 是边界层污染能否扩散的重要指示\n",
    "\n",
    "具体如何识别利用 TLogP 图 ，可以参考 [气象家园 TLogP 课件](http://bbs.06climate.com/forum.php?mod=viewthread&tid=28868&extra=&page=1 \"气象家园 TLogP 课件\")\n",
    "\n",
    "## 微波辐射计数据绘制 TLogP 图\n",
    "\n",
    "metpy 库基于 matplotlib 库提供了 TlogP 图绘制的接口 skew ，可以很方便地将高空数据绘制成 TLogP 图。\n",
    "\n",
    "该接口需要使用 气压作为 参数传入，微波辐射计采集的高空气象数据只有高度数据和地面气压数据，没有高空气压数据。\n",
    "\n",
    "可以参照气压高度递减率的方法 ，使用 metpy.calc.add_height_to_pressure 方法计算不同层的气压。\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a8d440b-9494-45ef-b0a0-5345275561a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 拆分微波辐射计数据，分为地面和高空两组数据\n",
    "sur_data = df[['DateTime', 'SurTem(℃)', 'SurHum(%)', 'SurPre(hPa)',\n",
    "               'Tir(℃)', 'Rain', 'CloudBase(km)', 'Vint(mm)', 'Lqint(mm)']]\n",
    "high_data = df.drop(['Record',  'SurTem(℃)', 'SurHum(%)', 'SurPre(hPa)',\n",
    "                     'Tir(℃)', 'Rain', 'CloudBase(km)', 'Vint(mm)', 'Lqint(mm)', 'QCflag'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201e1e2f-1d09-4eb7-a08f-1c66924f81d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import metpy.calc as mpcalc\n",
    "from metpy.units import units\n",
    "# 丢弃重复值\n",
    "sur_data = sur_data.drop_duplicates().reset_index(drop=True)\n",
    "# 计算地面露点温度\n",
    "sur_data['SurDwp'] = mpcalc.dewpoint_from_relative_humidity(\n",
    "    sur_data['SurTem(℃)'].values*units.degC,\n",
    "    sur_data['SurHum(%)'].values * units.percent).m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17db8fd8-7576-475e-96f7-cccf15a14cc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b544ccc2-ef1b-4a87-b9f0-f28e76504f2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data.columns  # 查看列名"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e6d70ec-8727-471c-94e8-01afc2e92b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 长表转横，把原先的不同廓线转换为同一行数据。\n",
    "high_data = high_data.pivot(index='DateTime', columns='dtype', values=high_data.columns[2:])\n",
    "\n",
    "# 并把高度由 km 转换成 m 单位，形成 _ 分隔的宽表\n",
    "high_data.columns = [al[1] +'_' +str(int(float(al[0].replace('(km)', '')) * 1000)) for al in high_data.columns]\n",
    "\n",
    "# 重置索引\n",
    "high_data = high_data.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96035e09-0be5-4ec8-9741-257eaf901340",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data  # 查看数据情况"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91ab6acb-d68e-49dc-b8ee-b3b0cd092077",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 宽表转长，这样就可以得到不同时间不同层的温度，为下一步计算和绘图做准备\n",
    "high_data = pd.wide_to_long(\n",
    "    high_data,\n",
    "    stubnames=[\n",
    "        \"RHU\",\n",
    "        \"TEM\",\n",
    "        \"WDen\",\n",
    "        \"WVDen\"],\n",
    "    i=\"DateTime\",\n",
    "    j=\"height\",\n",
    "    sep=\"_\").reset_index()\n",
    "high_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c34c94ad-389d-4d78-9a1c-7ecb50cd32e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data = pd.merge(\n",
    "    high_data, sur_data[['DateTime',  'SurTem(℃)', 'SurHum(%)', 'SurDwp', 'SurPre(hPa)']], on='DateTime', how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e0bdf2e-a96f-467a-a092-482432efea0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2a1980b-da51-48d6-a5bf-81c8a7ff62ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 计算露点温度 和 不同高度的气压\n",
    "# 使用 metpy.calc.dewpoint_from_relative_humidity(temperature, relative_humidity) 方法计算不同高度的露点温度\n",
    "# 使用 metpy.calc.add_height_to_pressure(pressure, height) 方法计算地面以上的气压\n",
    "\n",
    "high_data['DWP'] = mpcalc.dewpoint_from_relative_humidity(\n",
    "    high_data['TEM'].values*units.degC,\n",
    "    high_data['RHU'].values * units.percent).m\n",
    "high_data['PRS_HWC'] = mpcalc.add_height_to_pressure(\n",
    "    high_data['SurPre(hPa)'].values * units.hPa,\n",
    "    high_data['height'].values*units.meter).m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27cfa56d-693c-4a6e-a65c-86e742358d1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc0b0f16-0d73-4f79-89d1-7bdfa66015b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_data = high_data[high_data['DateTime'] == datetime(\n",
    "    2024, 12, 4, 8)].sort_values('height', ascending=True).reset_index(drop=True)\n",
    "show_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36dff348-dd69-43d2-b36f-bfcf10cf20de",
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy.plots import SkewT\n",
    "\n",
    "# 设置绘图参数\n",
    "fig = plt.figure(figsize=(18, 12))\n",
    "skew = SkewT(fig, rotation=0)#注意 rotation 表示图像旋转，默认为 30度 ，旋转表示斜温图\n",
    "\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats()\n",
    "skew.plot_mixing_lines()\n",
    "skew.ax.set_xlim(-80, 30)\n",
    "skew.ax.set_ylim(1000, 200)\n",
    "skew.ax.set_title('温度压力对数图', fontsize=18)\n",
    "\n",
    "# 微波辐射计的数据\n",
    "p = show_data['PRS_HWC'].values * units.hPa\n",
    "T = show_data['TEM'].values*units.degC\n",
    "Td = show_data['DWP'].values*units.degC\n",
    "\n",
    "# 绘制 温度层结曲线\n",
    "skew.plot(p, T, 'blue')\n",
    "# 绘制 露点层结曲线\n",
    "skew.plot(p, Td, 'green')\n",
    "\n",
    "\n",
    "# 绘制 LCL 抬升凝结高度\n",
    "lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])\n",
    "skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')\n",
    "\n",
    "# 计算气块上升的温度数据\n",
    "prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')\n",
    "# 绘制状态曲线\n",
    "skew.plot(p, prof, 'maroon', linewidth=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a6051ae-e352-4619-a7fe-f7b239c8ee66",
   "metadata": {},
   "source": [
    "## TLogP 图 叠加风廓线雷达的数据\n",
    "\n",
    "经过上图绘制，我们发现虽然已经绘制 三条探空数据曲线，但是与 TLogP 介绍的截图相比，还缺少风羽。\n",
    "\n",
    "虽然微波辐射计无法探测高空的风场，但是上一章节我们学习的 风廓线雷达 可以探测风场数据。\n",
    "\n",
    "往往我们把微波辐射计和风廓线雷达、及一些相关垂直遥感观测系统统称为 垂直观测系统。这些设备相互之间互相印证补充，起到了  $1+1>2$ 的效果。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320d1872-ee26-4099-b939-1acc3294ce88",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 由于 风廓线雷达采用 UTC 时间，因此采用00时的数据文件。\n",
    "fkxdi=WindProfileRadar.readSingleL3file(datapath+'WindProfileRadarFiles/20241204/WIND_ROBS/Z_RADA_I_77777_20241203200000_P_WPRD_LC_ROBS.TXT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d573ee3c-68fe-4a38-886e-788cb8b0f3e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fkx_df=fkxdi.Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c3bec6a-b928-48a9-aee3-fb2099710911",
   "metadata": {},
   "outputs": [],
   "source": [
    "fkx_df['SurPre(hPa)']=show_data['SurPre(hPa)'].values[0]#把地面气压加入，方便后续处理\n",
    "#计算各高度的气压\n",
    "fkx_df['PRS_HWC'] = mpcalc.add_height_to_pressure(\n",
    "    fkx_df['SurPre(hPa)'].values * units.hPa,\n",
    "    fkx_df['Sampling_Height'].values*units.meter\n",
    ").m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae5c73a6-fa1f-4f4b-ac57-f81599ec42cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fkx_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87ada22f-37c8-48df-94f7-3622fb952231",
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy.plots import SkewT\n",
    "# 设置绘图参数\n",
    "fig = plt.figure(figsize=(18, 12))\n",
    "skew = SkewT(fig, rotation=0)#注意 rotation 表示图像旋转，默认为 30度 ，旋转表示斜温图\n",
    "\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats()\n",
    "skew.plot_mixing_lines()\n",
    "skew.ax.set_xlim(-80, 30)\n",
    "skew.ax.set_ylim(1000, 200)\n",
    "skew.ax.set_title('温度压力对数图', fontsize=18)\n",
    "\n",
    "\n",
    "# 微波辐射计的数据\n",
    "p = show_data['PRS_HWC'].values * units.hPa\n",
    "T = show_data['TEM'].values*units.degC\n",
    "Td = show_data['DWP'].values*units.degC\n",
    "# 风廓线仪器的数据\n",
    "p2 = fkx_df['PRS_HWC'].values * units.hPa\n",
    "u = fkx_df['U_Wind_Speed'].values * units.knots\n",
    "v = fkx_df['V_Wind_Speed'].values * units.knots\n",
    "\n",
    "# 绘制 温度层结曲线\n",
    "skew.plot(p, T, 'blue')\n",
    "\n",
    "# 绘制 露点层结曲线\n",
    "skew.plot(p, Td, 'green')\n",
    "\n",
    "# 绘制 LCL 抬升凝结高度\n",
    "lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])\n",
    "skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')\n",
    "\n",
    "# 计算气块上升期的温度\n",
    "prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')\n",
    "# 绘制 状态曲线\n",
    "skew.plot(p, prof, 'maroon', linewidth=2)\n",
    "\n",
    "# 绘制 风羽\n",
    "skew.plot_barbs(p2, u, v);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc780f18-91fa-4bd8-b557-0ce063a1bf9f",
   "metadata": {},
   "source": [
    "叠加风羽的 TLogP 图更能分析高空冷暖平流输送，对灾害性天气的预测有更大的作用。\n",
    "\n",
    "TLogP 图也常常被称为 TLnP 图，在 气象 MICAPS 软件中又被称为 **温度压力对数图**。\n",
    "\n",
    "颜色为 blue 的曲线表示不同气压高度下的大气温度廓线，即**温度层结曲线**；\n",
    "\n",
    "颜色为 green  的曲线表示不同气压高度下 大气露点温度廓线，即**露点层结曲线**；\n",
    "\n",
    "颜色为 maroon  的曲线表示虚拟气块上升的温度状态曲线，即**状态曲线**；\n",
    "\n",
    "由于该日无强对流，因无法绘制 **CAPE （对流有效位能）** ；不过从图中我们可以发现 近地层 **露点层结曲线和温度层结曲线很接近**，说明该日**近地面相对湿度较高**，结合风羽看，**风力较小**，这样往往存在 **轻雾或大雾等天气**。\n",
    "\n",
    "实况也是如此，该日能见度自 07 时不断走低，最小能见度 200 米发生在 08 时至 09 时之间。\n",
    "\n",
    "当然有经验的气象预报员还可以从探空图中看出更多的信息，如果屏幕前的你是气象系统职工，可以上中国气象局远程培训网查找相应课程学习。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "528d7a4e-a17d-4373-8b27-e4d315ec0b2f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244e530f-8e0d-48d9-977c-bf5bd467279f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c752000-c45e-407a-904f-c90efc87463f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e4246a6-14a8-4d7c-b271-e83c18b7785e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e800d562-2e18-41b5-aed0-5245c1f52ca8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:qxfx3p11]",
   "language": "python",
   "name": "conda-env-qxfx3p11-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
